Load packages

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(RColorBrewer)

Variable importance graphic ~ Individual FRRs RF

SDI

2012

library(tidyverse)
# The only thing you should have to change in this script is the "path" 
# in list.files and the path of the variable match file. 
varmatch <- read.csv("./data/varmatch.csv")

tfiles <- list.files(path = "./results/",
                     pattern = "RF12_SDIimp_[[:digit:]]+.RDS", 
                     full.names = TRUE)

levels <- c("Mississippi Portal: 50.61%","Basin and Range: 34.86%","Fruitful Rim: 49.07%","Southern Seaboard: 48.50%","Eastern Uplands: 64.19%","Prairie Gateway: 51.71%", "Northern Great Plains: 54.81%", "Northern Crescent: 69.28%", "Heartland: 77.21%")

temp <- vector("list", length(tfiles))
for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::rename(IncPurity = tidyr::starts_with("IncPurity")) %>%
    dplyr::mutate(FRR = levels[i], 
                  Prop = IncPurity/sum(IncPurity)*100)
}

# Combine and match the custom variable names
rf12sdi_imp_frr2 <- data.table::rbindlist(temp) %>%
  dplyr::mutate(FRR = factor(FRR, levels = levels)) %>%
  dplyr::left_join(., varmatch, by = "varnames") %>%
  dplyr::mutate(category = factor(category, unique(varmatch$category))) %>%
  dplyr::arrange(category)

rf12sdi_imp_frr2$varfull <- factor(rf12sdi_imp_frr2$varfull, 
                                   levels = unique(rf12sdi_imp_frr2$varfull))

rf12sdi_imp_frr2 <- rf12sdi_imp_frr2 %>% group_by(FRR) %>% 
  mutate(stand = IncPurity/max(IncPurity)) %>% # standardized importance 
  arrange(desc(IncPurity))

# size indicates proportion of importance where most important = largest bubble and least important = smallest bubble
rf12sdi_imp_frrviz <- ggplot(rf12sdi_imp_frr2, aes(varfull,FRR)) +
  geom_point(aes(size=stand*8.4, colour = category), shape = 20,) +
  #scale_size(range = c(.1,15)) +
  scale_size_identity() +
  theme_minimal()+
  theme(axis.text.x=element_text(angle=55,hjust=1,vjust=1),axis.text.y=element_text(size=10),legend.text = element_text(size = 12),legend.title = element_blank(), legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5)))+
  scale_colour_brewer(palette = "Dark2") +
  ylab("") +
  xlab("") +
  ggtitle("SDI 2012")

rf12sdi_imp_frrviz

ggsave(filename = "./for_manuscript/rf12sdi_imp_frr.png", plot = rf12sdi_imp_frrviz, dpi = 300, width = 8.5, height = 6)

2017

varmatch <- read.csv("./data/varmatch.csv")

tfiles <- list.files(path = "./results/",
                     pattern = "RF17_SDIimp_[[:digit:]]+.RDS", 
                     full.names = TRUE)

levels <- c("Mississippi Portal: 9.58%*","Basin and Range: 33.88%","Fruitful Rim: 42.81%","Southern Seaboard: 61.72%","Eastern Uplands: 68.35%","Prairie Gateway: 47.93%", "Northern Great Plains: 58.64%", "Northern Crescent: 65.42%", "Heartland: 73.86%")

temp <- vector("list", length(tfiles))
for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::rename(IncPurity = tidyr::starts_with("IncPurity")) %>%
    dplyr::mutate(FRR = levels[i], 
                  Prop = IncPurity/sum(IncPurity)*100)
}

# Combine and match the custom variable names
rf17sdi_imp_frr2 <- data.table::rbindlist(temp) %>%
  dplyr::mutate(FRR = factor(FRR, levels = levels)) %>%
  dplyr::left_join(., varmatch, by = "varnames") %>%
  dplyr::mutate(category = factor(category, unique(varmatch$category))) %>%
  dplyr::arrange(category)

rf17sdi_imp_frr2$varfull <- factor(rf17sdi_imp_frr2$varfull, 
                                   levels = unique(rf17sdi_imp_frr2$varfull))

rf17sdi_imp_frr2 <- rf17sdi_imp_frr2 %>% group_by(FRR) %>% 
  mutate(stand = IncPurity/max(IncPurity)) %>% # standardized importance 
  arrange(desc(IncPurity))

# size indicates proportion of importance where most important = largest bubble and least important = smallest bubble
rf17sdi_imp_frrviz <- ggplot(rf17sdi_imp_frr2, aes(varfull,FRR)) +
  geom_point(aes(size=stand*8.4, colour = category), shape = 20,) +
  #scale_size(range = c(.1,15)) +
  scale_size_identity() +
  theme_minimal()+
  theme(axis.text.x=element_text(angle=55,hjust=1,vjust=1),axis.text.y=element_text(size=10),legend.text = element_text(size = 12),legend.title = element_blank(), legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5)))+
  scale_colour_brewer(palette = "Dark2") +
  ylab("") +
  xlab("") +
  ggtitle("SDI 2017")

rf17sdi_imp_frrviz

ggsave(filename = "./for_manuscript/rf17sdi_imp_frr.png", plot = rf17sdi_imp_frrviz, dpi = 300, width = 8.5, height = 6)

SIDI

2012

varmatch <- read.csv("./data/varmatch.csv")

tfiles <- list.files(path = "./results/",
                     pattern = "RF12_SIDIimp_[[:digit:]]+.RDS", 
                     full.names = TRUE)

levels <- c("Mississippi Portal: 49.82%","Basin and Range: 31.17%","Fruitful Rim: 35.02%","Southern Seaboard: 62.23%","Eastern Uplands: 67.01%","Prairie Gateway: 44.13%", "Northern Great Plains: 52.19%", "Northern Crescent: 68.20%", "Heartland: 70.19%")

temp <- vector("list", length(tfiles))
for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::rename(IncPurity = tidyr::starts_with("IncPurity")) %>%
    dplyr::mutate(FRR = levels[i], 
                  Prop = IncPurity/sum(IncPurity)*100)
}

# Combine and match the custom variable names
rf12sidi_imp_frr2 <- data.table::rbindlist(temp) %>%
  dplyr::mutate(FRR = factor(FRR, levels = levels)) %>%
  dplyr::left_join(., varmatch, by = "varnames") %>%
  dplyr::mutate(category = factor(category, unique(varmatch$category))) %>%
  dplyr::arrange(category)

rf12sidi_imp_frr2$varfull <- factor(rf12sidi_imp_frr2$varfull, 
                                   levels = unique(rf12sidi_imp_frr2$varfull))

rf12sidi_imp_frr2 <- rf12sidi_imp_frr2 %>% group_by(FRR) %>% 
  mutate(stand = IncPurity/max(IncPurity)) %>% # standardized importance 
  arrange(desc(IncPurity))

# size indicates proportion of importance where most important = largest bubble and least important = smallest bubble
rf12sidi_imp_frrviz <- ggplot(rf12sidi_imp_frr2, aes(varfull,FRR)) +
  geom_point(aes(size=stand*8.4, colour = category), shape = 20,) +
  #scale_size(range = c(.1,15)) +
  scale_size_identity() +
  theme_minimal()+
  theme(axis.text.x=element_text(angle=55,hjust=1,vjust=1),axis.text.y=element_text(size=10),legend.text = element_text(size = 12),legend.title = element_blank(), legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5)))+
  scale_colour_brewer(palette = "Dark2") +
  ylab("") +
  xlab("") +
  ggtitle("SIDI 2012")

rf12sidi_imp_frrviz

ggsave(filename = "./for_manuscript/rf12sidi_imp_frr.png", plot = rf12sidi_imp_frrviz, dpi = 300, width = 8.5, height = 6)

2017

varmatch <- read.csv("./data/varmatch.csv")

tfiles <- list.files(path = "./results/",
                     pattern = "RF17_SIDIimp_[[:digit:]]+.RDS", 
                     full.names = TRUE)

levels <- c("Mississippi Portal: 10.26%*","Basin and Range: 31.11%","Fruitful Rim: 39.76%","Southern Seaboard: 50.73%","Eastern Uplands: 65.29%","Prairie Gateway: 48.24%", "Northern Great Plains: 56.29%", "Northern Crescent: 68.20%", "Heartland: 70.19%")

temp <- vector("list", length(tfiles))
for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::rename(IncPurity = tidyr::starts_with("IncPurity")) %>%
    dplyr::mutate(FRR = levels[i], 
                  Prop = IncPurity/sum(IncPurity)*100)
}

# Combine and match the custom variable names
rf17sidi_imp_frr2 <- data.table::rbindlist(temp) %>%
  dplyr::mutate(FRR = factor(FRR, levels = levels)) %>%
  dplyr::left_join(., varmatch, by = "varnames") %>%
  dplyr::mutate(category = factor(category, unique(varmatch$category))) %>%
  dplyr::arrange(category)

rf17sidi_imp_frr2$varfull <- factor(rf17sidi_imp_frr2$varfull, 
                                   levels = unique(rf17sidi_imp_frr2$varfull))

rf17sidi_imp_frr2 <- rf17sidi_imp_frr2 %>% group_by(FRR) %>% 
  mutate(stand = IncPurity/max(IncPurity)) %>% # standardized importance 
  arrange(desc(IncPurity))

# size indicates proportion of importance where most important = largest bubble and least important = smallest bubble
rf17sidi_imp_frrviz <- ggplot(rf17sidi_imp_frr2, aes(varfull,FRR)) +
  geom_point(aes(size=stand*8.4, colour = category), shape = 20,) +
  #scale_size(range = c(.1,15)) +
  scale_size_identity() +
  theme_minimal()+
  theme(axis.text.x=element_text(angle=55,hjust=1,vjust=1),axis.text.y=element_text(size=10),legend.text = element_text(size = 12),legend.title = element_blank(), legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5)))+
  scale_colour_brewer(palette = "Dark2") +
  ylab("") +
  xlab("") +
  ggtitle("SIDI 2017")

rf17sidi_imp_frrviz

ggsave(filename = "./for_manuscript/rf17sidi_imp_frr.png", plot = rf17sidi_imp_frrviz, dpi = 300, width = 8.5, height = 6)

SIDI TRANSF

2012 TRANSF

varmatch <- read.csv("./data/varmatch.csv")

tfiles <- list.files(path = "./results/",
                     pattern = "RF12_SIDITRANSFimp_[[:digit:]]+.RDS", 
                     full.names = TRUE)

levels <- c("Mississippi Portal: 36.78%","Basin and Range: 33.93%","Fruitful Rim: 42.33%","Southern Seaboard: 52.14%","Eastern Uplands: 58.67%","Prairie Gateway: 48.43%", "Northern Great Plains: 54.37%", "Northern Crescent: 66.96%", "Heartland: 71.48%")

temp <- vector("list", length(tfiles))
for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::rename(IncPurity = tidyr::starts_with("IncPurity")) %>%
    dplyr::mutate(FRR = levels[i], 
                  Prop = IncPurity/sum(IncPurity)*100)
}

# Combine and match the custom variable names
rf12siditransf_imp_frr2 <- data.table::rbindlist(temp) %>%
  dplyr::mutate(FRR = factor(FRR, levels = levels)) %>%
  dplyr::left_join(., varmatch, by = "varnames") %>%
  dplyr::mutate(category = factor(category, unique(varmatch$category))) %>%
  dplyr::arrange(category)

rf12siditransf_imp_frr2$varfull <- factor(rf12siditransf_imp_frr2$varfull, 
                                   levels = unique(rf12siditransf_imp_frr2$varfull))

rf12siditransf_imp_frr2 <- rf12siditransf_imp_frr2 %>% group_by(FRR) %>% 
  mutate(stand = IncPurity/max(IncPurity)) %>% # standardized importance 
  arrange(desc(IncPurity))

# size indicates proportion of importance where most important = largest bubble and least important = smallest bubble
rf12siditransf_imp_frrviz <- ggplot(rf12siditransf_imp_frr2, aes(varfull,FRR)) +
  geom_point(aes(size=stand*8.4, colour = category), shape = 20,) +
  #scale_size(range = c(.1,15)) +
  scale_size_identity() +
  theme_minimal()+
  theme(axis.text.x=element_text(angle=55,hjust=1,vjust=1),axis.text.y=element_text(size=10),legend.text = element_text(size = 12),legend.title = element_blank(), legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5)))+
  scale_colour_brewer(palette = "Dark2") +
  ylab("") +
  xlab("") +
  ggtitle("SIDI 2012")

rf12siditransf_imp_frrviz

ggsave(filename = "./for_manuscript/rf12siditransf_imp_frr.png", plot = rf12siditransf_imp_frrviz, dpi = 300, width = 8.5, height = 6)

2017 TRANSF

varmatch <- read.csv("./data/varmatch.csv")

tfiles <- list.files(path = "./results/",
                     pattern = "RF17_SIDITRANSFimp_[[:digit:]]+.RDS", 
                     full.names = TRUE)

levels <- c("Mississippi Portal: 13.13%*","Basin and Range: 32.77%","Fruitful Rim: 31.30%","Southern Seaboard: 60.46%","Eastern Uplands: 64.96%","Prairie Gateway: 47.81%", "Northern Great Plains: 51.33%", "Northern Crescent: 60.70%", "Heartland: 65.34%")

temp <- vector("list", length(tfiles))
for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::rename(IncPurity = tidyr::starts_with("IncPurity")) %>%
    dplyr::mutate(FRR = levels[i], 
                  Prop = IncPurity/sum(IncPurity)*100)
}

# Combine and match the custom variable names
rf17siditransf_imp_frr2 <- data.table::rbindlist(temp) %>%
  dplyr::mutate(FRR = factor(FRR, levels = levels)) %>%
  dplyr::left_join(., varmatch, by = "varnames") %>%
  dplyr::mutate(category = factor(category, unique(varmatch$category))) %>%
  dplyr::arrange(category)

rf17siditransf_imp_frr2$varfull <- factor(rf17siditransf_imp_frr2$varfull, 
                                   levels = unique(rf17siditransf_imp_frr2$varfull))

rf17siditransf_imp_frr2 <- rf17siditransf_imp_frr2 %>% group_by(FRR) %>% 
  mutate(stand = IncPurity/max(IncPurity)) %>% # standardized importance 
  arrange(desc(IncPurity))

# size indicates proportion of importance where most important = largest bubble and least important = smallest bubble
rf17siditransf_imp_frrviz <- ggplot(rf17siditransf_imp_frr2, aes(varfull,FRR)) +
  geom_point(aes(size=stand*8.4, colour = category), shape = 20,) +
  #scale_size(range = c(.1,15)) +
  scale_size_identity() +
  theme_minimal()+
  theme(axis.text.x=element_text(angle=55,hjust=1,vjust=1),axis.text.y=element_text(size=10),legend.text = element_text(size = 12),legend.title = element_blank(), legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5)))+
  scale_colour_brewer(palette = "Dark2") +
  ylab("") +
  xlab("") +
  ggtitle("SIDI 2017")

rf17siditransf_imp_frrviz

ggsave(filename = "./for_manuscript/rf17siditransf_imp_frr.png", plot = rf17siditransf_imp_frrviz, dpi = 300, width = 8.5, height = 6)

RICH

2012

varmatch <- read.csv("./data/varmatch.csv")

tfiles <- list.files(path = "./results/",
                     pattern = "RF12_RICHimp_[[:digit:]]+.RDS", 
                     full.names = TRUE)

levels <- c("Mississippi Portal: 54.54%","Basin and Range: 49.86%","Fruitful Rim: 68.70%","Southern Seaboard: 60.96%","Eastern Uplands: 56.10%","Prairie Gateway: 38.67%", "Northern Great Plains: 48.39%", "Northern Crescent: 68.55%", "Heartland: 46.05%")

temp <- vector("list", length(tfiles))
for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::rename(IncPurity = tidyr::starts_with("IncPurity")) %>%
    dplyr::mutate(FRR = levels[i], 
                  Prop = IncPurity/sum(IncPurity)*100)
}

# Combine and match the custom variable names
rf12rich_imp_frr2 <- data.table::rbindlist(temp) %>%
  dplyr::mutate(FRR = factor(FRR, levels = levels)) %>%
  dplyr::left_join(., varmatch, by = "varnames") %>%
  dplyr::mutate(category = factor(category, unique(varmatch$category))) %>%
  dplyr::arrange(category)

rf12rich_imp_frr2$varfull <- factor(rf12rich_imp_frr2$varfull, 
                                   levels = unique(rf12rich_imp_frr2$varfull))

rf12rich_imp_frr2 <- rf12rich_imp_frr2 %>% group_by(FRR) %>% 
  mutate(stand = IncPurity/max(IncPurity)) %>% # standardized importance 
  arrange(desc(IncPurity))

# size indicates proportion of importance where most important = largest bubble and least important = smallest bubble
rf12rich_imp_frrviz <- ggplot(rf12rich_imp_frr2, aes(varfull,FRR)) +
  geom_point(aes(size=stand*8.4, colour = category), shape = 20,) +
  #scale_size(range = c(.1,15)) +
  scale_size_identity() +
  theme_minimal()+
  theme(axis.text.x=element_text(angle=55,hjust=1,vjust=1),axis.text.y=element_text(size=10),legend.text = element_text(size = 12),legend.title = element_blank(), legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5)))+
  scale_colour_brewer(palette = "Dark2") +
  ylab("") +
  xlab("") +
  ggtitle("RICH 2012")

rf12rich_imp_frrviz

ggsave(filename = "./for_manuscript/rf12rich_imp_frr.png", plot = rf12rich_imp_frrviz, dpi = 300, width = 8.5, height = 6)

2017

varmatch <- read.csv("./data/varmatch.csv")

tfiles <- list.files(path = "./results/",
                     pattern = "RF17_RICHimp_[[:digit:]]+.RDS", 
                     full.names = TRUE)

levels <- c("Mississippi Portal: 37.66%","Basin and Range: 29.91%","Fruitful Rim: 71.32%","Southern Seaboard: 61.86%","Eastern Uplands: 50.80%","Prairie Gateway: 48.11%", "Northern Great Plains: 31.09%", "Northern Crescent: 61.93%", "Heartland: 55.7%")

temp <- vector("list", length(tfiles))
for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::rename(IncPurity = tidyr::starts_with("IncPurity")) %>%
    dplyr::mutate(FRR = levels[i], 
                  Prop = IncPurity/sum(IncPurity)*100)
}

# Combine and match the custom variable names
rf17rich_imp_frr2 <- data.table::rbindlist(temp) %>%
  dplyr::mutate(FRR = factor(FRR, levels = levels)) %>%
  dplyr::left_join(., varmatch, by = "varnames") %>%
  dplyr::mutate(category = factor(category, unique(varmatch$category))) %>%
  dplyr::arrange(category)

rf17rich_imp_frr2$varfull <- factor(rf17rich_imp_frr2$varfull, 
                                   levels = unique(rf17rich_imp_frr2$varfull))

rf17rich_imp_frr2 <- rf17rich_imp_frr2 %>% group_by(FRR) %>% 
  mutate(stand = IncPurity/max(IncPurity)) %>% # standardized importance 
  arrange(desc(IncPurity))

# size indicates proportion of importance where most important = largest bubble and least important = smallest bubble
rf17rich_imp_frrviz <- ggplot(rf17rich_imp_frr2, aes(varfull,FRR)) +
  geom_point(aes(size=stand*8.4, colour = category), shape = 20,) +
  #scale_size(range = c(.1,15)) +
  scale_size_identity() +
  theme_minimal()+
  theme(axis.text.x=element_text(angle=55,hjust=1,vjust=1),axis.text.y=element_text(size=10),legend.text = element_text(size = 12),legend.title = element_blank(), legend.position = "bottom") +
  guides(color = guide_legend(override.aes = list(size = 5)))+
  scale_colour_brewer(palette = "Dark2") +
  ylab("") +
  xlab("") +
  ggtitle("RICH 2017")

rf17rich_imp_frrviz

ggsave(filename = "./for_manuscript/rf17rich_imp_frr.png", plot = rf17rich_imp_frrviz, dpi = 300, width = 8.5, height = 6)

PDP of the same variables across all regions

Check distribution of data

load("./data/div_response.RData") 

hist(div_17$chem) #log transform

hist(div_17$fert) #log transform

hist(div_17$perc_cl)

hist(div_17$BV4)

hist(div_17$BV15)

hist(div_17$gvt_prog) #log transform

hist(div_17$acres_per_op)

hist(div_17$T_OC)

SDI 2017

Chem

library(tidyverse)
library(cowplot)
library(sp)

load("./data/div_response.RData")

orig <- div_17
county <- readRDS("./data/county.RDS")
levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

tfiles <- list.files(path = "./results/",
                     pattern = "chem_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

pdp_chem <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

# Log Scale Version
tbreaks = c(0, 10, 20, 40, 80, 160, 320, 640)
tbreaks2 = c(0, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024)
ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)
#going to be different for each variable 
#look at lower part and read it regularly/normally 
#at the upper end, what happens when you use a lot of chemical
#use log scale to actually interpret what is happening 

chem_pp_sdi2017 <- ggplot(pdp_chem) +   
  geom_line(aes(x = log2(x + 1), y = y, color = VAR), size = 1.2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) +
  scale_y_continuous(breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Chemical input (2017)*") + 
  theme_classic()+ 
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

# chem_pp_sdi2017 +   geom_rug(orig, mapping = aes(x = chem, color = FRR_NAME)) 

orig$FRR_NAME <- factor(orig$FRR_NAME, levels = levels)

chem_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = log2(chem + 1), y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) + 
  ylab("") + xlab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic()+  
  theme(legend.position = "none")


jpeg("./for_manuscript/RFpdp_SDI2017_chem.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(chem_pp_sdi2017, chem_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2

Fert

levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

tfiles <- list.files(path = "./results/",
                     pattern = "fert_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

pdp_fert <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

# Log Scale Version
tbreaks = c(0, 10, 20, 40, 80, 160, 320, 640)
tbreaks2 = c(0, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024)
ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)

#going to be different for each variable 
#look at lower part and read it regularly/normally 
#at the upper end, what happens when you use a lot of chemical
#use log scale to actually interpret what is happening 

fert_pp_sdi2017 <- ggplot(pdp_fert) +   
  geom_line(aes(x = log2(x + 1), y = y, color = VAR), size = 1.2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) +
  scale_y_continuous(limits = c(0.6, 1.3), breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Fertilizer input (2017)*") + 
  theme_classic() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

# chem_pp_sdi2017 +   geom_rug(orig, mapping = aes(x = chem, color = FRR_NAME)) 

orig$FRR_NAME <- factor(orig$FRR_NAME, levels = levels)
fert_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = log2(fert + 1), y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) + 
  ylab("") + xlab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic() + 
  theme(legend.position = "none")


jpeg("./for_manuscript/RFpdp_SDI2017_fert.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(fert_pp_sdi2017, fert_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2

Government Programs

levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

summary(orig$gvt_prog)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.00903  3.26836  8.68279 11.57501 17.18242 82.95022
tfiles <- list.files(path = "./results/",
                     pattern = "gp_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

pdp_gp <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

# Log Scale Version
tbreaks = c(0, 5, 10, 20, 40, 60, 120, 340)
tbreaks2 = c(0, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024)
ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)
#going to be different for each variable 
#look at lower part and read it regularly/normally 
#at the upper end, what happens when you use a lot of chemical
#use log scale to actually interpret what is happening 

gp_pp_sdi2017 <- ggplot(pdp_gp) +   
  geom_line(aes(x = log2(x + 1), y = y, color = VAR), size = 1.2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) +
  scale_y_continuous(limits = c(0.6, 1.3), breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Government programs (2017)*") + 
  theme_classic() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

# chem_pp_sdi2017 +   geom_rug(orig, mapping = aes(x = chem, color = FRR_NAME)) 

orig$FRR_NAME <- factor(orig$FRR_NAME, levels = unique(pdp_gp$VAR))
gp_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = log2(gvt_prog + 1), y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) + 
  ylab("") + xlab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic() + 
  theme(legend.position = "none")


jpeg("./for_manuscript/RFpdp_SDI2017_gp.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(gp_pp_sdi2017, gp_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2

Precipitation seasonality

levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

summary(orig$BV15) #coefficient of variation
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   35.51   51.06   61.22   66.41   80.28  149.66
tfiles <- list.files(path = "./results/",
                     pattern = "ps_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

#Log-Scale Version
tbreaks2 = c(32, 64, 128, 256)
ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)

pdp_ps <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

ps_pp_sdi2017 <- ggplot(pdp_ps) +   
  geom_line(aes(x = log2(x + 1), y = y, color = VAR), size = 1.2) + 
scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2)+
  scale_y_continuous(limits = c(0.6, 1.3), breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Precipitation seasonality (2017)*") + 
  theme_classic() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

orig$FRR_NAME <- factor(orig$FRR_NAME, levels = levels)

ps_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = log2(BV15 + 1), y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) +
  ylab("") + xlab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2)+
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic() + 
  theme(legend.position = "none")

jpeg("./for_manuscript/RFpdp_SDI2017_ps.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(ps_pp_sdi2017, ps_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2

Temperature seasonality

levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

summary(orig$BV4) #SD*100
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   317.8   796.4   902.5   902.1  1009.7  1311.2
tfiles <- list.files(path = "./results/",
                     pattern = "ts_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)

pdp_ts <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

ts_pp_sdi2017 <- ggplot(pdp_ts) +   
  geom_line(aes(x = x, y = y, color = VAR), size = 1.2) + 
scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", limits = c(300, 1350)) +
  scale_y_continuous(limits = c(0.6, 1.3), breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Temperature seasonality (2017)") + 
  theme_classic() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

# chem_pp_sdi2017 +   geom_rug(orig, mapping = aes(x = chem, color = FRR_NAME)) 

orig$FRR_NAME <- factor(orig$FRR_NAME, levels = levels)
ts_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = BV4, y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) + 
  ylab("") + xlab("") + 
  scale_x_continuous("", limits = c(300, 1350)) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic() + 
  theme(legend.position = "none")

jpeg("./for_manuscript/RFpdp_SDI2017_ts.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(ts_pp_sdi2017, ts_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2

Percent cropland

levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

tfiles <- list.files(path = "./results/",
                     pattern = "cl_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)

pdp_cl <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

cl_pp_sdi2017 <- ggplot(pdp_cl) +   
  geom_line(aes(x = x, y = y, color = VAR), size = 1.2) + 
scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", limits = c(0, 100)) +
  scale_y_continuous(limits = c(0.6, 1.3), breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Percent cropland (2017)") + 
  theme_classic() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

# chem_pp_sdi2017 +   geom_rug(orig, mapping = aes(x = chem, color = FRR_NAME)) 

orig$FRR_NAME <- factor(orig$FRR_NAME, levels = levels)

cl_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = perc_cl, y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) + 
  ylab("") + xlab("") + 
  scale_x_continuous("", limits = c(0, 100)) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic() + 
  theme(legend.position = "none")

jpeg("./for_manuscript/RFpdp_SDI2017_cl.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(cl_pp_sdi2017, cl_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2

Topsoil organic carbon

summary(orig$T_OC)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3484  0.9403  1.0663  1.3766  1.4461 34.8717
levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

tfiles <- list.files(path = "./results/",
                     pattern = "toc_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

tbreaks2 = c(0, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024)
ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)
#going to be different for each variable 
#look at lower part and read it regularly/normally 
#at the upper end, what happens when you use a lot of chemical
#use log scale to actually interpret what is happening 

pdp_toc <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

toc_pp_sdi2017 <- ggplot(pdp_toc) +   
  geom_line(aes(x = log2(x + 1), y = y, color = VAR), size = 1.2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) +
  scale_y_continuous(limits = c(0.6, 1.3), breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Topsoil organic carbon (2017)*") + 
  theme_classic() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))


orig$FRR_NAME <- factor(orig$FRR_NAME, levels = levels)
toc_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = log2(T_OC + 1), y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) + 
  ylab("") + xlab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic() + 
  theme(legend.position = "none")


jpeg("./for_manuscript/RFpdp_SDI2017_toc.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(toc_pp_sdi2017,toc_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2

Farm size

summary(orig$acres_per_op)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     3.0    50.0    84.0   170.6   150.0  4550.0
levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

tfiles <- list.files(path = "./results/",
                     pattern = "fs_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

tbreaks2 = c(0, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024)
ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)
#going to be different for each variable 
#look at lower part and read it regularly/normally 
#at the upper end, what happens when you use a lot of chemical
#use log scale to actually interpret what is happening 

pdp_fs <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

fs_pp_sdi2017 <- ggplot(pdp_fs) +   
  geom_line(aes(x = log2(x + 1), y = y, color = VAR), size = 1.2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) +
  scale_y_continuous(limits = c(0.6, 1.3), breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Farm size (2017)*") + 
  theme_classic() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

# chem_pp_sdi2017 +   geom_rug(orig, mapping = aes(x = chem, color = FRR_NAME)) 

orig$FRR_NAME <- factor(orig$FRR_NAME, levels = levels)
fs_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = log2(acres_per_op + 1), y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) + 
  ylab("") + xlab("") + 
  scale_x_continuous("", breaks = log2(tbreaks2 + 1), labels = tbreaks2) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic() + 
  theme(legend.position = "none")


jpeg("./for_manuscript/RFpdp_SDI2017_fs.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(fs_pp_sdi2017, fs_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2

Percent pastureland (excluding cropland)

levels = c("Mississippi Portal", "Basin and Range","Fruitful Rim", "Southern Seaboard","Eastern Uplands","Prairie Gateway","Northern Great Plains","Northern Crescent", "Heartland")

levels_rev = rev(levels)

tfiles <- list.files(path = "./results/",
                     pattern = "pe_sdi17_[[:digit:]]+.RDS", full.names = TRUE)

temp <- vector("list", length(tfiles))

for(i in 1:length(tfiles)){
  temp[[i]] <- readRDS(tfiles[i]) %>%
    dplyr::mutate(VAR = levels_rev[i])
}

ybreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3)

pdp_pe <- data.table::rbindlist(temp) %>%
  dplyr::mutate(VAR = factor(VAR, levels = levels))

pe_pp_sdi2017 <- ggplot(pdp_pe) +   
  geom_line(aes(x = x, y = y, color = VAR), size = 1.2) + 
scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  ylab("") + 
  scale_x_continuous("", limits = c(0, 100)) +
  scale_y_continuous(limits = c(0.6, 1.3), breaks = ybreaks, labels = ybreaks)+
  ggtitle("SDI - Percent pastureland (2017)") + 
  theme_classic() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

# chem_pp_sdi2017 +   geom_rug(orig, mapping = aes(x = chem, color = FRR_NAME)) 

orig$FRR_NAME <- factor(orig$FRR_NAME, levels = levels)

pe_bar_2017 <- ggplot(orig) + 
  geom_boxplot(aes(x = perc_pe, y = FRR_NAME, color = FRR_NAME), outlier.shape = "|", outlier.stroke = 3) + 
  ylab("") + xlab("") + 
  scale_x_continuous("", limits = c(0, 100)) + 
  scale_colour_manual(values=c("#CCCC99", "#336633", "#99CC99", "#646667","#CCCCCC", "#4A7DA5","#99CCFF","#FFC300","#CC6633")) +
  theme_classic() + 
  theme(legend.position = "none")


jpeg("./for_manuscript/RFpdp_SDI2017_pe.jpg",
    width = 6,
    height = 5,
    res = 300,
    units = "in")
plot_grid(pe_pp_sdi2017, pe_bar_2017, ncol=1, align="v", rel_heights = c(2, 1))
dev.off()
## quartz_off_screen 
##                 2